How Much is a Coffee Worth to Me?

I moved to Melbourne last year, which, along with being a questionable career choice, turned an already troubling dependence on caffeine into a behavioural disorder.

I bought a Nespresso machine—telling myself I would save money on coffee—and it has been a huge success. Now, with the help of my Nespresso machine I have more than doubled my coffee consumption. The coffee pod coffee doesn’t count, you see. Yes I’ll come down for a coffee, I haven’t had mine yet.

This leads me to wonder how much of my family’s future I am eroding to bean-dust, every time I head downstairs for another medium latte, thanks—the absent droning sound I make every morning between the hours of 9 and 10 am. Each milky-upper is costing me $6. So assuming I work 250 days in a year… let’s see now:

COFFEE_COST = 6.0
WORKING_DAYS = 250

annual_cost = COFFEE_COST * WORKING_DAYS

It’s $1500.0. You can check the code if you want to make sure I got that.

I’m not devastated by that number—it’s not a small amount, but it’s not enough to pay for much in Australia, either. Maybe a less masochistic person would leave the matter there, but it begs the question of what I might be doing with that money if I wasn’t spending it all on morning browns to straighten out my morning frowns.

Opportunity Cost, and Futures Lost…

To estimate at how much this behaviour is costing me, I’ve put together a model that compares two strategies:

1. Nick buys a coffee every day
2. Nick invests the price of a coffee every day into the market, instead

After 20 years, we’ll see how Nick the milk dribbler is doing, compared to Nick the slightly more financially responsible investor.

Basic Model

Before things get more complicated, let’s have a look at how these strategies perform assuming an 8% annual return on the market:

def create_date_array(start_date: datetime, num_months: int):
    start = np.datetime64(start_date.replace(day=1), 'M')
    month_array = start + np.arange(num_months)
    return month_array.astype("datetime64[D]")

def annual_to_monthly_rate(annual_rate: float) -> float:
    return (1 + annual_rate) ** (1 / 12) - 1

def calculate_compound_returns(
    investments: np.ndarray, 
    annual_rate: float
) -> np.ndarray:
    num_months = len(investments)
    monthly_return = annual_to_monthly_rate(annual_rate)
    returns = np.zeros(num_months)

    for i in range(num_months):
        months_invested = np.arange(num_months - i)
        returns[i:] += investments[i] * (1 + monthly_return) ** months_invested

    return returns


YEARS = 20
ANNUAL_RATE = 0.08
DAYS_IN_MONTH = 30.437      # Accounting for gap years

num_months = YEARS * 12
months = np.arange(num_months)
dates = create_date_array(datetime.today(), num_months=num_months)

investing = np.full_like(dates, COFFEE_COST * DAYS_IN_MONTH).astype(float)
spending = -1 * investing
returns = calculate_compound_returns(
    investments=investing,
    annual_rate=ANNUAL_RATE
)

df = pd.DataFrame({
    "Date": dates,
    "Coffee": np.cumsum(spending),
    "Investing": np.cumsum(investing),
    "Returns": returns - np.cumsum(investing),
}).round(2)
flat_df = df.melt(
    id_vars=["Date"], 
    var_name="Case", 
    value_name="$"
)

def plot_balance_area_chart(df: pd.DataFrame, title: str):
    return (
        alt.Chart(df)
        .mark_area(opacity=0.75)
        .encode(
            alt.X("Date:T", title=""),
            alt.Y("$:Q"),
            alt.Color(
                "Case:N",
                legend=alt.Legend(orient="top-left")
            ),
            tooltip=[
                alt.Tooltip("Date:T", title=""),
                alt.Tooltip("$:Q", title="$", format="$,.2f"),
                alt.Tooltip("Case:N", title="")
            ],
            order=alt.Order("Case:N")
        )
        .properties(
            title=title,
            height=PLOT_HEIGHT,
            width=PLOT_WIDTH
        )
    )

chart = plot_balance_area_chart(
    flat_df, title="Comparison of Investing vs Coffee Drinking"
)
chart.show()

total_gains = int(df.iloc[-1]["Investing"] + df.iloc[-1]["Returns"])
total_losses = int(df.iloc[-1]["Coffee"])

Unless I’ve screwed up the maths here, things are already looking very grim. Nick the dribbler is down $43680 in 20 years, whereas Nick the investor is giving his kid a very nice $103557 graduation bonus.

Devaluing the Future

Ok, but what about inflation, and the time value of money, and all that finance stuff? asks Nick the dribbler, pretending he isn’t an idiot.

Well the time value of money is a bit of a flimsy concept when the present money is being spent on hot milk that will cool down in minutes and expire in a matter of days, but let’s see.

Assuming 2.5% inflation annually, and a discount rate of 13.5%, which look like this:

def create_exponential_curve(
    months: np.ndarray, 
    annual_rate: float
) -> np.ndarray:
    monthly_rate = annual_to_monthly_rate(annual_rate)
    return np.exp(months * monthly_rate)


DISCOUNT_RATE = -0.135
INFLATION_RATE = 0.025

discount_curve = create_exponential_curve(months, DISCOUNT_RATE)
inflation_curve = create_exponential_curve(months, INFLATION_RATE)

def plot_factor_curves(discount_curve, inflation_curve) -> alt.Chart:
    return (
        alt.Chart(
            pd.DataFrame({
                "Date": dates,
                "Inflation": inflation_curve,
                "Discount": discount_curve,
                "Net Adjustment": inflation_curve * discount_curve
            })
            .melt(
                id_vars=["Date"], 
                var_name="Curve",
                value_name="Factor"
            )
        )
        .mark_line()
        .encode(
            alt.X("Date:T"),
            alt.Y("Factor:Q"),
            alt.Color(
                "Curve:N",
                legend=alt.Legend(orient="top-left"),
            ),
            tooltip=[
                alt.Tooltip("Date:T", title=""),
                alt.Tooltip("Factor:Q", format=".3f"),
                alt.Tooltip("Curve:N", title="")
            ],
        )
        .properties(
            height=PLOT_HEIGHT,
            width=PLOT_WIDTH
        )
    )

chart = plot_factor_curves(discount_curve, inflation_curve)
chart.show()

… Investor Nick’s returns look like this:

discount_df = pd.DataFrame({
    "Date": dates,
    "Coffee": np.cumsum(spending) * inflation_curve * discount_curve,
    "Investing": np.cumsum(investing) * inflation_curve * discount_curve,
    "Returns": (returns - np.cumsum(investing)) * inflation_curve * discount_curve,
}).round(2)

flat_discount_df = discount_df.melt(
    id_vars=["Date"], 
    var_name="Case", 
    value_name="$"
)

chart = plot_balance_area_chart(
    flat_discount_df, 
    title="(Real, Discounted at 13.5%): Comparison of Investing vs Coffee Drinking"
)
chart.show()

…So what does that mean?

If present me values future dollars at a discount rate of 13.5%, I am basically claiming that money on the table today can be converted into a 13.5% return in a year, through my… savvy investing strategy.

That means that future dollars are worth less to me (not worthless, worth less). The sooner I have those dollars, the sooner I can put them to work to make me that return.

And compounding this out to 20 years from now implies that a dollar then is worth about 5% of its value to me today.

This applies to dollars spent as well as dollars invested—future expenses are smaller when discounted too.

So what discount rate could justify Nick the dribbler’s behaviour? Well, he isn’t investing that money, which means he must really, really value that 3 minute mouth experience, which would imply a very aggressive discount rate, probably in excess of 100%.

discount_curve = create_exponential_curve(months, annual_rate=-1)

chart = plot_factor_curves(discount_curve, inflation_curve)
chart.show()

discount_df = pd.DataFrame({
    "Date": dates,
    "Coffee": np.cumsum(spending) * inflation_curve * discount_curve,
    "Investing": np.cumsum(investing) * inflation_curve * discount_curve,
    "Returns": (returns - np.cumsum(investing)) * inflation_curve * discount_curve,
}).round(2)

flat_discount_df = discount_df.melt(
    id_vars=["Date"], 
    var_name="Case", 
    value_name="$"
)

chart = plot_balance_area_chart(
    flat_discount_df, 
    title="(Real, Discounted at 100%): Comparison of Investing vs Coffee Drinking"
)
chart.show()

So, any money beyond a year time horizon is worth jack shit to Nick the dribbler. He lives by the froth, dies by the froth, neglects his financial future by the froth.

From a maximize money standpoint this isn’t rational, and realistically speaking I would expect an individual to have a discount rate of a little over some risk-free return—like treasury bonds at around 4.5%—up to around 9%, depending on their risk tolerance. Here’s 6%:

discount_curve = create_exponential_curve(months, annual_rate=-0.06)
discount_df = pd.DataFrame({
    "Date": dates,
    "Coffee": np.cumsum(spending) * inflation_curve * discount_curve,
    "Investing": np.cumsum(investing) * inflation_curve * discount_curve,
    "Returns": (returns - np.cumsum(investing)) * inflation_curve * discount_curve,
}).round(2)

flat_discount_df = discount_df.melt(
    id_vars=["Date"], 
    var_name="Case", 
    value_name="$"
)

chart = plot_balance_area_chart(
    flat_discount_df, 
    title="(Real, Discounted at 6%): Comparison of Investing vs Coffee Drinking"
)
chart.show()

Stochastic Model

But you can’t just assume a 8.0% year on year return like that! Market performance isn’t guaranteed! protests Nick the dribbler, milk froth and steam spilling from his maw as if Smaug the dragon had recently left the Lonely Mountain, put on some weight and taken up residence in South Melbourne.

And fair enough! From a risk perspective it’s worth examining what range of outcomes someone is exposed to when thinking about saving for the future. It’s also not as daunting a task as it might first appear, especially if we settle for a relatively simple to implement price model.

First, we need data.

Historical Returns

Here is the last 15 years of ASX200 market data:

import yfinance as yf

price_df = (
    yf.Ticker("^AXJO")
    .history(period="15y")
    .reset_index()
)

chart = (
    alt.Chart(price_df)
    .mark_line(color="red", strokeWidth=1)
    .encode(
        alt.X("Date:T"),
        alt.Y(
            "Close:Q",
            scale=alt.Scale(zero=False)
        )
    )
    .properties(
        title='ASX 200 Closing Price - 15 Years',
        height=PLOT_HEIGHT,
        width=PLOT_WIDTH,
    )
)

chart.show()

Since the model invests into “the market”, I’ll just take the ASX200 as representative of market movements to keep things simple.

For our forecast, we’ll need to transform this data into something a bit more convenient, the log returns, defined as:

\[ \begin{align} \ln R_t = \ln \frac{P_t}{P_{t-1}} \end{align} \]

Where \(P_t\) is the closing price on a given day, and \(P_{t-1}\) is the closing price on the previous day.

price_df["Returns"] = (
    np.log(
        price_df["Close"] / price_df["Close"].shift(1)
    )
)

returns_chart = (
    alt.Chart(price_df)
    .mark_line(color="limegreen", strokeWidth=1)
    .encode(
        alt.X("Date:T"),
        alt.Y("Returns:Q")
    )
    .properties(
        title='ASX 200 Log Returns',
        height=PLOT_HEIGHT,
        width=PLOT_WIDTH,
    )
)

returns_chart.show()

Because we are (more or less) taking the derivative of price—the degree to which prices change day to day—we end up with a nice, relatively stationary signal of the market. Removing the time axis gives us a dataset that I can fit a distribution to:

params = stats.norm.fit(price_df["Returns"].dropna())
dist = stats.norm(*params)
x_range = np.linspace(
    price_df["Returns"].min(), 
    price_df["Returns"].max(),
    250
)

chart = stat_plots.hist_dist_plot(
    price_df,
    col="Returns",
    scipy_dist=dist,
    title="Normal Distribution Fitted to Returns",
    hist_color="limegreen"
)

chart.show()

The normal distribution used above does a pretty poor job of capturing the data. A Student’s T distribution should better capture the tails:

params = stats.t.fit(price_df["Returns"].dropna())
dist = stats.t(*params)
x_range = np.linspace(
    price_df["Returns"].min(), 
    price_df["Returns"].max(),
    250
)

chart = stat_plots.hist_dist_plot(
    price_df,
    col="Returns",
    scipy_dist=dist,
    title="Student's T Distribution Fitted to Returns",
    hist_color="limegreen"
)

chart.show()

Much better.

At this point it looks like we have a good-enough fit, but it is worth plotting samples back on a time-axis so we can directly compare what our model produces to actual historical returns:

sample_df = pd.DataFrame({
    "Date": price_df["Date"],
    "Simulated Returns": dist.rvs(size=len(price_df["Date"]))
})

sample_chart = (
    alt.Chart(sample_df)
    .mark_line(color="orange", strokeWidth=1)
    .encode(
        alt.X("Date:T"),
        alt.Y("Simulated Returns:Q")
    )
)

combined_chart = (sample_chart + returns_chart).resolve_axis(x="shared", y="shared")
combined_chart.show()

… and it does okay. The elephant in the room is COVID-19. Independent samples from a distribution can’t capture clustered periods of volatility like that. The samples are also a bit fuzzier in general—there is more volatility in the simulated market than was seen in reality. Overall I give the model a:

This could be my career catch phrase.

Geometric Brownian Motion

Now that we have a distribution, we can simulate forecasts from it.

First, since the distribution is fitted to daily returns and our forecast will be monthly, an adjustment is needed:

def convert_t_dist_to_monthly(t_dist, days_per_month=21):
    df, mu, sigma = t_dist.args[0], t_dist.args[1], t_dist.args[2]

    # Monthly transformation
    mu_monthly = days_per_month * mu
    sigma_monthly = sigma * np.sqrt(days_per_month)

    return stats.t(df=df, loc=mu_monthly, scale=sigma_monthly)

monthly_dist = convert_t_dist_to_monthly(dist)

Right, now the fun stuff.

It turns out that you can cook up a stochastic forecast called a Geometric Brownian Motion (GBM) model quite easily with a numpy magic spell:

sims = 1000
initial_balance = COFFEE_COST

def simulate_monthly_gbm(
    monthly_returns_dist,
    initial_balance: float,
    num_months: int,
    sims: int = 100
) -> np.ndarray:

    return initial_balance * (
        np.exp(
            np.cumsum(
                monthly_returns_dist.rvs(
                    size=(num_months, sims)
                ),
                axis=0
            )
        )
    )

balance = simulate_monthly_gbm(
    monthly_dist,
    initial_balance=COFFEE_COST,
    num_months=num_months,
    sims=sims
)

In words, this is the exponentiated cumulative sum of log returns, or excumsumlogret.

“EXCUMSUUUMMLOGREEEEET”

chart = finance_plots.simulated_portfolio_plot(
    simulated_balance=balance,
    dates=dates,
    title="Coffee Simulations",
    sims=100
)

chart.show()

When you hold out your coffee cup and say the magic word, future worlds spout forth, in which the $6.0 purchase of that coffee is instead invested into the market.

Looking at the end point of all those universes 20 years from now, we get a distribution of that coffee’s final value:

df = pd.DataFrame({
    "Sim": np.arange(sims),
    "$": balance[-1, :]
})

chart = (
    alt.Chart(df)
    .mark_bar()
    .encode(
        alt.X("$:Q", bin=alt.Bin(maxbins=250)),
        alt.Y(
            "count()",
            title="",
            axis=alt.Axis(labels=False, ticks=False, title=None),
        )
    )
    .properties(
        title='Final Coffee Value',
        height=PLOT_HEIGHT,
        width=PLOT_WIDTH,
    )
)

chart.show()

percentiles = [1, 10, 50, 90, 99]

df = pd.DataFrame(
    {
       "Coffee Value": np.percentile(df["$"], q=percentiles),
    }, 
    index=[f"P{100 - p}" for p in percentiles]
).round()

df.T

	P99	P90	P50	P10	P1
Coffee Value	23.0	45.0	107.0	249.0	519.0

Simulated Investments

Now that we have a way of simulating forecasts, I can run the investing strategy (putting the price of a cup of coffee into the market each day) through the model to see how the portfolio might perform over the 20 year period.

def simulate_compounded_returns(
    investments: np.ndarray,
    monthly_returns_dist,
    num_months: int,
) -> np.ndarray:

    returns = np.exp(
        monthly_returns_dist.rvs(size=(num_months, sims))
    )

    portfolio_values = np.zeros_like(returns)
    for i in range(num_months):
        # Initial investment
        if i == 0:
            portfolio_values[i] = investments[i] * returns[i]

        # New investment + previously invested
        else:
            portfolio_values[i] = (
                (portfolio_values[i-1] + investments[i]) * returns[i]
            )

    return portfolio_values

portfolio = simulate_compounded_returns(
    investments=investing,
    monthly_returns_dist=monthly_dist,
    num_months=num_months
)

chart = finance_plots.simulated_portfolio_plot(
    portfolio,
    dates,
    title="Simulated Investment Portfolio",
    sims=100
)

chart.show()

df = pd.DataFrame({
    "Sim": np.arange(sims),
    "$": portfolio[-1, :]
})

chart = (
    alt.Chart(df)
    .mark_bar()
    .encode(
        alt.X("$:Q", bin=alt.Bin(maxbins=250)),
        alt.Y(
            "count()",
            title="",
            axis=alt.Axis(labels=False, ticks=False, title=None),
        )
    )
    .properties(
        title='Final Portfolio Value',
        height=PLOT_HEIGHT,
        width=PLOT_WIDTH,
    )
)

chart.show()

And if we discount at the conservative rate of 6%:

discounted_portfolio = portfolio * discount_curve[:, np.newaxis]

chart = finance_plots.simulated_portfolio_plot(
    discounted_portfolio,
    dates,
    title="(Discounted) Simulated Investment Portfolio",
    sims=100
)
chart.show()

df = pd.DataFrame({
    "Sim": np.arange(sims),
    "$": discounted_portfolio[-1, :]
})

chart = (
    alt.Chart(df)
    .mark_bar()
    .encode(
        alt.X("$:Q", bin=alt.Bin(maxbins=250)),
        alt.Y(
            "count()",
            title="",
            axis=alt.Axis(labels=False, ticks=False, title=None),
        )
    )
    .properties(
        title='(Discounted) Final Portfolio Value',
        height=PLOT_HEIGHT,
        width=PLOT_WIDTH,
    )
)

chart.show()

percentiles = [1, 10, 50, 90, 99]

df = pd.DataFrame(
    {
       "Nominal Returns": np.percentile(portfolio[-1, :], q=percentiles),
       "Discounted Returns": np.percentile(discounted_portfolio[-1, :], q=percentiles)
    }, 
    index=[f"P{100 - p}" for p in percentiles]
).round()

df.T

	P99	P90	P50	P10	P1
Nominal Returns	82134.0	134836.0	268181.0	502817.0	950964.0
Discounted Returns	24027.0	39444.0	78452.0	147090.0	278188.0

So we can’t know for sure how the market will perform over the next 20 years, but it looks like a good bet to drink less coffee and invest more. There’s quite a lot of positive skew in the model, so there is a lot of upside potential in investing that money.

… That feels like a very underwhelming conclusion, after all that effort.

What’s the point of this, again?

The Good

Even though all models are wrong and every prediction is contingent on a set of assumptions and all that… I think the power of doing this is to stoke the imagination about what might be possible beyond the typical point estimate like 20.54% return over 10 years, especially when it comes to thinking about exposure—to risks as well as opportunities.

Having a single prediction doesn’t help you to think about the range of things that might happen, but a stochastic model like the simple forecast above does. Now you can think about how likely it is you’ll lose all your money, or how likely it is that curbing your coffee habit will make you a millionaire, or anything in between.

And what about the implications of that range of outcomes? The lower tail predictions of the model might have extremely dire consequences that lead you not to act, for fear of the repercussions. The upper tail might be so attractive that it becomes obvious that the exposure is worth the entry price, even if the “expectation” doesn’t look that great. This kind of reasoning isn’t possible without uncertainty quantification.

The robustness of the model is always going to be a concern, but it seems to me a lot less flimsy than a single prediction. I’m not sure quantitative predictions about the future make much sense at all unless they quantify uncertainty somehow. Why draw a particular scenario out of a hat when there are an infinite set of others to choose from?

The Bad

But this also makes stochastic models way more dangerous—having a distribution as output might lead you or your colleagues to believe you’ve actually captured the possibilities objectively. You haven’t. What you have is the formalized consequences of your assumptions—in this case, a simplistic representation of past behaviour of the market propagated into the future—and that can still be useful, but it’s not the same as knowing what will actually happen. COVID19 is not in your model. World War III is not in your model.

Point estimates don’t carry that pretense. They’re dumb, but at least they’re honest about it. Or rather we are more honest about them.

This is a real tension I’ve felt in my career. The whole point of providing estimates is to help people make better decisions, but the uncertainty in those estimates is not true—it’s the output of a model. Models aren’t reality. This becomes obvious when you change a model hyperparameter and see all the estimates change—like an epistemology witch crawling out of the screen and whispering in your ear, it’s just a model, stupid. But that might not be as obvious when viewed from the outside.

“But cross validation is so bori—AAAAAHHHH!”

So when I hand off my estimates… are we all just pretending that we “know” the uncertainty now? If the managers understand they should take the estimates with a “grain of salt”… how many grains of salt? At some point human intuition steps in and takes over again, which diminishes the point of the exercise.

The Coffee

Take this coffee example. The model says that on paper Nick the investor is a lot better off, only having lost the value of the pleasure of a coffee… and the relationship building that goes along with it… and the chance conversations that could have really impacted his career… every single morning… for 20 years…

Maybe the real investments were the friends we made along the way…

And there it is! My intuitions (biased as they are) are wrestling with the implications of the model. I know it is not very sophisticated, I know that life is about more than maximizing financial gain, and I know that I can’t quantify those things, so it is inevitably neglecting many factors I care about.

So when faced with that doubt, my decision to drink coffee loops back to aesthetics, or how I want to live my life, not a quantitative, value maximizing rational one.

So modelling things is useful, but decision making is a lot more nuanced and human than numerical modelling and decision theory would like to make it.

… I probably should cut down, though.

Thanks for reading.